Domain Growth, Budding, and Fission in Phase Separating Self- Assembled Fluid 

Bilayers 
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A systematic investigation of the phase separation dynamics in self-assembled multi-component 
bilayer fluid vesicles and open membranes is presented. We use large-scale dissipative particle 
dynamics to explicitly account for solvent, thereby allowing for numerical investigation of the effects 
of hydrodynamics and area-to-volume constraints. In the case of asymmetric lipid composition, we 
observed regimes corresponding to coalescence of flat patches, budding, vesiculation and coalescence 
of caps. The area-to- volume constraint and hydrodynamics have a strong influence on these regimes 
and the crossovers between them. In the case of symmetric mixtures, irrespective of the area-to- 
volume ratio, we observed a growth regime with an exponent of 1/2. The same exponent is also 
found in the case of open membranes with symmetric composition. 

PACS numbers: 87.16.-b, 64.75. +g, 68.05.Cf 



I. INTRODUCTION 

Biomembranes are fascinating self- assembled quasi- 
two-dimensional complex fluids composed essentially of 
phospholipids and cholesterol. The primary roles of 
biomembranes are the separation between the inner and 
outer environments of the cell or inner organelles, and the 
support of an amazing specialized protein-based machin- 
ery which is crucial for a variety of physiological func- 
tions, transmembrane transport and structural integrity 
of the cell Jj. Many recent experiments demonstrated 
that biomembranes of eukariotic cells are laterally orga- 
nized into small nanoscopic domains, called rafts, which 
are rich in sphingomyelin and cholesterol. The higher 
content of cholesterol in rafts is due to the fact that 
the acyl chains in sphingomyelin are mainly saturated, 
thereby promoting their interaction with cholesterol. Al- 
though, it is largely believed that this inplane organi- 
zation is essential for a variety of physiological func- 
tions such as signaling, recruitment of specific proteins 
and endocytosis elucidation of the fundamental is- 
sues including the mechanisms leading to the formation 
of lipid rafts, their stability, and finite size remain elusive. 
Clearly, raft formation in biomembranes is complicated 
by the presence of many non-equilibrium mechanisms. In 
view of this, it is important to understand the equilibrium 
phase behavior and the kinetics of fluid multicomponent 
lipid membranes before attempts are made to find the ef- 
fects of more complex mechanisms that may be involved 
in the formation and stability of lipid rafts. 

The dynamics of in-plane demixing in multicomponent 
lipid membranes is richer than their counterparts in Eu- 
clidean three- or two-dimensional systems. This is largely 
due to (i) the strong coupling between the lipid composi- 
tion and the membrane conformation, (ii) the difference 
between the viscosities of the lipid bilayer and that of the 
embedding fluid, and (iii) the area-to- volume constraint, 
maintained by a gradient in osmotic pressure across the 



membrane. As a result, various growth regimes, may 
be observed in multicomponent lipid membranes. Phase 
separation dynamics in multicomponent vesicles follow- 
ing a quench from a single homogeneous phase to the 
two-phase liquid-liquid coexistence region of the phase 
diagram has previously been considered by means of a 
generalized time-dependent Ginzburg-Landau model on 
a non-Euclidean closed surface 0, 13- Limitations im- 
posed by the parameterization of surface deformations 
did not allow for budding in this study. More recent 
simulations using a dynamic triangulation Monte Carlo 
model |H0il3' predicted dynamics which are much more 
complex than that in Euclidean surfaces. In particular, 
these simulations showed that for symmetric composition 
of binary mixtures, at intermediate times, the dynamics 
is characterized by the presence of the familiar labyrinth- 
like spinodal pattern. At later times, in the presence 
of curvature composition coupling, these patterns break 
up leading to isolated islands. At still later times, and 
in the case of tensionless membranes, these domains re- 
shape into buds connected to the parent membrane by 
very narrow necks. Further domain growth proceeds via 
the Brownian diffusion of these buds, and their coales- 
cence. It is important to remark that both the gener- 
alized time-dependent Ginzburg-Landau model and the 
Monte Carlo dynamic triangulation model (1) do not ac- 
count for the solvent, and are therefore purely dissipative, 
(2) cannot account for the constraint on the volume en- 
closed by the vesicle, (3) conserve the topology of the 
membrane throughout the simulation, and therefore do 
not account for fission or fusion processes. A model that 
accounts of these effects is clearly warranted in order to 
compare with experiments. In order to achieve this we 
used a model based on the dissipative particle dynamics 
(DPD) approach 0. 

On the experimental side, recent studies using ad- 
vanced techniques such as two-photon fluorescence and 
confocal microscopy, performed on giant unilamellar vesi- 
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cles (GUVs) composed of dioleoylphosphatidylcholine, 
sphingomyelin, and cholesterol, showed the existence of 
lipid domains 0. More recent experiments on ternary 
mixtures composed of saturated and unsaturated phos- 
phatidylcholine and cholesterol 0, 0, ^| saw the ex- 
istence of liquid-liquid coexistence over a wide range of 
compositions, an indication that liquid-liquid coexistence 
in lipid membranes is ubiquitous to a wide variety of 
ternary lipid mixtures. It must be noted, however, that 
domains observed in these experiments are comparable 
to the size of the vesicle (micron-scale) , orders of magni- 
tudes larger than rafts in biomcmbranes. Some of these 
experiments reported structures with many, more-or-less, 
curved domains. But these are more akin to caps than 
fully developed buds. An important question that arises 
out of these studies is: (1) What role does the solvent 
hydrodynamics and the volume constraints play during 
coarsening in multicomponent fluid vesicles? (2) Does 
the topology change drastically alter the kinetic pathway 
predicted by simulations that does not allow for topology 
changes? 



In this paper, we present results from extensive nu- 
merical simulations of self-assembled lipid vesicles and 
open membranes using the dissipative particle dynamics 
approach. We specifically investigated the phase sepa- 
ration dynamics in lipid membranes following a quench 
from the one phase region to the fluid-fluid coexistence 
region of the phase diagram. The lipid membrane is 
composed of self-assembled lipid particles in an explicit 
solvent, thus accounting for hydrodynamic effects. Fur- 
thermore, the parameters of the model are such that the 
membrane is impermeable to the solvent, thus allowing 
us to investigate the effect of area-to- volume ratio on the 
dynamics. We specifically investigated the effects of (1) 
area-to- volume ratio, (2) line tension, and (3) lipids com- 
position on the dynamics. We find that the path, through 
which dynamics proceeds, depends on the area-to- volume 
ratio and composition. In off-critical quenches, in partic- 
ular, the dynamics proceeds via the coalescence of small 
flat patches at intermediate times, followed by their bud- 
ding and vesiculation. At late times, domain growth pro- 
ceeds via coalescence of caps remaining on the vesicle. 
Crossovers between these regimes are strongly affected 
by the area-to-volume ratio and line tension. In the case 
of critical quenches, domain growth proceeds via dynam- 
ics similar to that in Euclidian two-dimensional fluids. 
That is, in this case the effect on the embedding fluid on 
the coarsening process seems to be not very obvious. We 
check this by also by simulating critical quenches in open 
membranes with different projected areas. 



This article is organized as follows: in Sec. II, the 
model and simulation technique are presented. In Sec. 
Ill, results of our simulations are presented. Finally, we 
summarize and conclude in Sec. IV. 
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FIG. 1: Snapshots of a phase separating vesicle with low ex- 
cess area, and (I)b ~ 0.5 at (a) t = 100 r, (b) 500 r, (c) 1000 r, 
(d) 2000 r, (e) 3000 t, and (f) 4000 t. 



II. MODEL AND METHOD 

In this work, the dynamics of phase de-mixing of a 
two-component lipid mixture in an explicit solvent is in- 
vestigated using dissipative particle dynamics (DPD) ap- 
proach. DPD was first introduced by Hoogerbrugge and 
Koelman more than a decade ago and was cast in 
its present form about five years later [T^ Il5l |. DPD 
is reminiscent of molecular dynamics, but is more ap- 
propriate for the investigation of generic properties of 
macromolecular systems. The use of soft repulsive in- 
teractions in DPD, allows for larger integration time in- 
crements than in a typical molecular dynamics simula- 
tion using Lennard- Jones interactions. Thus time and 
length scales much larger than those in atomistic molec- 
ular dynamics simulations can be probed by the DPD 
approach. Furthermore, DPD uses pairwise random and 
dissipative forces between neighboring particles, which 
are interrelated through the ffuctuation-dissipation the- 
orem. The pairwise nature of these forces ensures local 
conservation of momentum, a necessary condition for cor- 
rect long-range hydrodynamics [l^ . 

The system is composed of simple solvent particles (de- 
noted as w), and two types of complex lipid particles 
(A and B). A lipid particle is modeled as a flexible am- 
phiphilic chain with one hydrophilic particle, mimicking a 
lipid head group, and three hydrophobic particles, mim- 
icking the lipid acyl-group. More complicated lipid struc- 
tures and artifacts due to the choice of simulation param- 
eters have been investigated by other groups 0, ^| ■ 
These details are not expected to affect the qualitative 
nature of the results reported here. 
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The heads of the A and B hpids are denoted by Ha 
and Hb, and their respective tails are denoted by and 
tB- For simphcity, we focus in this study on the case 
where the interactions are symmetric under the exchange 
between A and B hpids. Thus this model does not con- 
tain any explicit coupling between curvature and compo- 
sition. The time evolution of the position and velocity 
of each dpd particle i, denoted by (ri,Vi), are governed 
by Hamilton's equations of motion. The three pairwise 
forces are given by 

Fzf^ = a^juj(r,j)n.,j, (1) 

F^f^ = 7^i'^^(''^J)(^■^J • Vjj)n,j, (2) 

F|f = a,,(At)i/2^(r,,)0,,n,„ (3) 

where = rj - r^, tiij = r^ /ry , and = - v^. % 
is a symmetric random variable satisfying 

= 0, (4) 

{O^J{t)0kl{t')) = [5,k5,i+5a5,k)5{t-t'). (5) 

with i ^ j and k ^ I. In Eq.©, At is the iteration time 
step. The weight factor is chosen as 

= I for r > re 

where Tc is the interactions cutoff radius. The choice of 
Lo in Eq. @ ensures that the interactions are all soft and 
repulsive. The equations of motion of particle i are given 
by 

dr,{t) = w,{t)dt, (7) 

^ 3 3 

+$:F(f(do^/^), (8) 

i 

where m is the mass of a single DPD particle. Here, 
for simplicity, masses of all types of dpd particles are 
supposed to be equal. Assuming that the system is in 
a heat bath at a temperature T, the parameters and 
aij in Eqs. ((21 and |3J) are related to each other by the 
fluctuation-dissipation theorem, 

7,, = 4-/2fcBr. (9) 

The parameters, , of the conservative forces are specif- 
ically chosen as 

/ hA tA w hB tB \ 

hA 25 200 25 oab 200 

tA 200 25 200 200 a^s 

~ w 25 200 25 25 200 ' ^ ' 

he ttAB 200 25 25 200 

V tB 200 aAB 200 200 25 / 
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FIG. 2: Snapshots of a phase separating vesicle in case-II. 
The time sequence of the snapshots is the same as in Fig. 1. 

where e and the cutoff radius, r^, set the energy and 
length scales, respectively. The effect of line tension is 
studied by varying the parameter aAB- The integrity of 
a lipid particle is ensured via an additional simple har- 
monic interaction, between consecutive particles, whose 
force is given by 

FMli = -C(l-^M+i/&)nM+i, (11) 

where we set, for the spring constant and the preferred 
bond length, values C — lOOe and b — 0.45rc, respec- 
tively. 

In our simulations, we used cr = 3.0 {e^m/riy^^. Most 
of the simulations were performed at fc^T = e, and a fluid 
density p — 3.Qt\7^. The iteration time was chosen to be 
At = 0.05t (2g, with the time scale r = {mr'l/ef/'^. 
The equations of motion are integrated using velocity- 
Verlet algorithm |23| . The total number of lipid particles 
used was 16 000, and both cases of closed vesicles and 
open membranes were simulated. In the case of a closed 
vesicle, the box size is (80 x 80 x 80)rj? corresponding to 
a total number of 1 536 000 dpd-particle. In the case of 
open membranes, we consider box sizes of Lx = Ly > Lz, 
with Lz = 40rc, such that the fluid density is equal to 
that in the case of closed vesicles. Open membranes with 
different tensions are simulated by varying , such that 
the number of lipid particles and fluid density are kept 
constant. Note that = 40rc is still much larger than 




FIG. 3: Net interface length as a function of time in case-I 
(blue curve) and case II (red curve). The green curve corre- 
sponds to case II with low line tension (see text for details of 
corresponding interactions). The dashed lines and the dot- 
dashed line have slopes of 0.3 and 4/9, respectively. In the 
inset, the number of vesiculated domains is shown as a func- 
tion of time for case-II with high line tension (corresponding 
to the red curve) . The arrows point to the time regime during 




FIG. 4: Number of clusters on the vesicle as a function of 
time. Curve colors correspond to those in Fig. |3] The slope 
of the dashed line is 2/3. Note that in the case of high excess 
area, clusters that have vesiculated are excluded. 



the thickness of the bilayer. Periodic boundary condi- 
tions are applied in all directions for both cases of closed 
and open membranes. 

Previous simulations based on DPD models have 
shown that open bilayers and closed vesicles can be self- 
assembled US 113' These studies, however, were 



FIG. 5: Mean square of the distance of the A-lipids from the 
center of mass. The top and bottom curves correspond to 
V = 0.17 and 0.14, respectively. 



performed on smaller systems than those in the present 
study. In order to save computer time, we prepare our 
vesicles starting from an almost closed configuration, 
composed from a single type of lipid. This approach al- 
lows for the equilibration of the lipid surface coverage in 
both leaflets through the diffusion of the lipids via the 
rim of the open vesicle. We find that the vesicle typically 
closes within about 50 000 DPD time steps. Once the 
vesicle is closed, we also found that within our param- 
eters, the number of solvent particles inside the vesicle 
and the numbers of lipid particles in the inner and outer 
leaflets remain constant throughout the simulation. A 
vesicle composed of 16 000 lipid chains, prepared as indi- 
cated above, is found to be nearly spherical and contain 
about 138 500 solvent particles inside it. Vesicles with ex- 
cess area (high area-to-volume ratio) are then prepared 
by transferring solvent particles from the core of the vesi- 
cles, prepared as indicated above, to the outer region, 
such that the fluid density is kept constant. An open 
membrane is prepared by placing a bilayer parallel to 
the a;y-plane at z = Lz/'i. The membrane is let to equi- 
librate until fluctuations of its height attain saturation. 
After equilibration of the closed vesicle or the open mem- 
brane, the phase separation process is triggered through 
an instantaneous change of a fraction of the A-lipids to 
B-lipids such that their composition is equal to (pB- This 
mimics a quench from a homogeneous state to the two- 
phase region. 

We have performed systematic simulations in which 
following parameters were varied: 

(i) The strength of the repulsive interaction between A 
and B lipids, aAB, hi order to infer the effect of line ten- 
sion, A, between A and B domains. This parameter was 
varied for both open and closed membranes. 

(ii) The compositions of the two lipids, (j)A and (pB = 



1— We considered the cases of (I>a = 0.5 and 0.3. This 
parameter was varied for both closed and open mem- 
branes. 

(iii) The area-to-volume ratio, i^, in the case of closed 
vesicles, defined here a.s v = [Nhead + Ntaii) /N^, where 
Nhead, Ntaii are the total numbers of head and tail dpd 
particles, respectively, and N^, is the total number of 
solvent particles inside the vesicle. In the case of open 
membranes, the projected area, effectively plays the role 
of area-to-volume ratio in closed vesicle. 

In the presentation and discussion of our results, we 
use the following labels to indicate the parameters used 
for different simulated systems 



system 


V 


aAB 




^(100) 


0.462 


100 


0.3 




0.567 


100 


0.3 




0.567 


86 


0.3 


_4(68) 


0.567 


68 


0.3 


/|(50) 


0.567 


50 


0.3 


^(100) 


0.462 


100 


0.5 


^(100) 


0.567 


50 


0.5 


^1 


0.567 


50 


0.5 



The interaction parameters of the present model are se- 
lected such that the membrane is impermeable to solvent 
particles. This implies that the number of solvent par- 
ticles inside closed vesicles is constant, thereby allowing 
us to effectively investigate the effect of area-to- volume 
ratio. In experiments, this parameter is controlled via 
the osmotic pressure. 

III. RESULTS AND DISCUSSION 

For the parameter values mentioned above, our model 
membrane does not show any flip-flop motion of the 
lipids, i.e. throughout the simulation time the number 
of lipids in the upper and lower membrane are the same. 
We also find that the coupling between the compositions 
of the two leaflets is found to be very strong. This is 
not surprising, considering the fact that we have chosen 

A. Case of closed vesicles with = 0.3 

In Fig. 1, snapshots of closed vesicles configurations in 
the case of 0b = 0.3 and with small area-to- volume ratio, 
corresponding to system are shown. This figure 

shows that the phase separation process after a quench 
to the two-phase region proceeds in a manner similar to 
that in Euclidian systems, i.e., through the formation of 
small domains, and their coalescence in time. During 
the early stages of the dynamics, domains have average 
curvatures that are equal to the surrounding majority 
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FIG. 6: Snapshot sequences for closed vesicles with excess 
area parameter u = 0.17. (a), (b) and (c) correspond to 
aAB = lOOe, 68e and 50e, respectively. Snapshots from top 
to bottom correspond to t = lOOr, lOOOr, 2000r, and 4000r, 
respectively. 



component, an indication that during the early stages 
of the phase separation process, the curvature is decou- 
pled from the composition. This is due to the fact that 
the compositions in the two leaflets are equal, and the 
A- and em B-lipids have identical architectures, leading 
to a decoupling between curvature and composition. As 
time evolves, the interface tension starts to assert and 
the domains curve very slightly from the majority com- 
ponent, while the vesicle becomes more spherical, an im- 
plication that an increase in tension has occurred. In 
Fig. 2, snapshots corresponding to (j)B ~ 0.3, but with a 
high excess area parameter , ^2^"*^^ are shown. The high 
excess area is obtained by removing 25610 solvent parti- 
cles from the core of the vesicle and putting them in the 
outer region, such that the fluid density is maintained. 
Notice, here again, that during the early stages, domains 
have same curvature as the majority component, imply- 
ing the decoupling between curvature and composition 
during these stages. At about 400r, B-domains start to 
curve away from the average vesicle's curvature. 

Domain growth is monitored through both the net in- 
terface length, L{t), and the average domain size, R{t), 
as calculated from the cluster size distribution. In order 
to see how the interfacial length can be used as a measure 




FIG. 7: Net interface length versus time for the case of closed 
vesicles with excess area parameter, v = 0.17, for different 
values of line tension. Curves from top to bottom correspond 
to aAB ~ 50e, 68e and lOOe, respectively. 



FIG. 8: The mean square of the distance of the A-lipid heads 
from the center of mass, (Ar)^ versus time. Curves from top 

and 



to bottom correspond to systems .Aj^"', J^^^'' 



(100) 



,(100) 
■^2 1 



of domain coarsening, consider a closed vesicle formed of 
f^it^ B-domains at some time with an average linear 
size i?. Furthermore, let us consider the case where the 
domains are circular. The net interface length is then 
given by Lit) = 27rN{t)R{t). Since the net amount of 
the B-component, is conserved, we also have for the area 
occupied by the B-component, Ab = it N (t) R'^ (t) . We 
therefore have 

L{t) ^ AB/R{t), (12) 
and the number of B-domains is therefore given by 

N{t)^ As/R^it). (13) 

In Fig. 3, L{t) versus time is shown for the case of (pA = 

0. 3 and for low and high area-to-volume parameter. We 
notice from this figure that, during intermediate times, 

1. e. t < 400t, the interfacial length is independent of the 
area-to- volume ratio, and has the form 

L(t)^t-", (14) 

with the growth exponent, a « 0.3. The number of do- 
mains as calculated from the cluster size distribution is 
shown in Fig. 4, which again shows that the number of 
clusters is independent of the area-to- volume ratio at in- 
termediate times, and that 

N{t) ^ (15) 

with /3 « 2q; « 2/3, in agreement with Eq. (|13l) . 

A growth exponent, a = 1/3, in Euclidian multicom- 
ponent systems is usually attributed to the evaporation- 
condensation mechanism as explained by the classical 
theory of Lifshitz and Slyozov (23|. In a fluid system 



such as ours, domain growth can also be the result of 
the motion of the entire domains themselves resulting in 
collision between domains leading to coalescence. 

Two domains coalesce if they travel a distance, l{t) de- 
termined by the average area on the membrane occupied 
by a domain. This is given by 

l{tf - A/N{t). (16) 

Moreover, assuming that the domains perform a Brown- 
ian walk before their collision, the collision time should 
obey l{t)'^ ^ Dfjt. In the case of an isolated two- 
dimensional fluid, as can be calculated from the Stoke's 
formula for the drag on a circular domain due to the sur- 
rounding fluid, the diffusion coefficient is independent of 
the domain size. However, in our case, the drag expe- 
rienced by the domains results mainly from the three- 
dimensional embedding fluid, leading to Da ^ 1/R. Us- 
ing this fact, we have the time dependence of domain size 
as 

i?(t) - (17) 
and the number of clusters 

in good agreement with our numerical results. 

Irrespective of the initial excess area, the very late 
times configuration is always that of a tensed vesicle. 
This can be seen from the mean square of of the po- 
sitions of the A-lipid heads from the center of mass of 
the vesicle, (Ar)^. This is shown in Fig. 5 for the case 

of GAB = lOOe, for systems and A^2"^\ A large 



m 




FIG. 9: The number of vesiculated buds versus time for sys- 
tern 4'°°^ 



(Ar)^ implies a floppy vesicle with a lot of excess area. 
Notice that (Ar)^ decreases rapidly after about 400r, 
eventually reaching a value equal to that without excess 
area. The fast decrease in (Ar)^ is due to the reshaping 
of the domains into caps, the budding of some of them 
and their vesiculation. The buds, once formed, are found 
to vesiculate within a short period of time, of order lOr. 
We confirm this by performing simulations of a single 
B-domain occupying 12% of the total area of a vesicle 
with excess area. The fission mechanism of a bud itself 
is a very interesting phenomena, and has recently been 
investigated in detail using DPD simulations [2^ . 

As a result of vesiculation, the parent vesicle looses 
most of its excess area, leading it to acquire a more spher- 
ical shape, as shown in Fig. 2, and implied by Fig. 5. This 
induced lateral tension prevents the domains from fur- 
ther capping. Further domain coarsening may therefore 
proceed via the coalescence of these caps. 



B. Effect of line tension on the phase separation in 
closed vesicles with (j>B = 0.3 



The effect of line tension on the dynamics of domain 
growth in closed vesicles with (pB = 0.3 is investigated by 
performing simulations at different values of gab = lOOe, 
86e, 68e, and 50e in the presence of large excess area. 

In Fig. 6, snapshots for the cases of gab = lOOe, 
68e and 50e are shown for comparison. These snapshots 
clearly show that the line tension play an important role 
on the dynamics. Corresponding interfacial lengths as a 
function of time are shown in Fig. 7. This figure again 
shows that budding is delayed as the line tension is in- 
creased. We also notice that while in the system with 
O'AB = lOOe, budding of domains is followed by their 
vesiculation, in the systems with uab — 86e and 68e, 





FIG. 10: Snapshot sequences of closed vesicles with (j>B = 0.5. 



(I), (II) and (III) correspond to systems C{ 



(100) ^(100) 



and 



€2'^'^^ , respectively. Snapshots from top to bottom correspond 
to t = lOOr, 200r, 400t and 700r, respectively. 



very few buds vesiculate. No vesiculation was observed 
in the case of low line tension {qab = 50e). 

In order to gain further insight into the effect of line 
tension on domains caping, let us consider a circular do- 
main of B-phase with area a, surrounded by a sea of 
A-phase, on a tensionless membrane. Let c be the abso- 
lute value of the mean curvature of this domain which is 
assumed to be uniform. The free energy of the domain 
is therefore given by 



£a = 2Kac -|- XI, 



(19) 



where k and A are the bending modulus and line tension, 
respectively, and I is the perimeter of the domains, given 
by 



I = 27r 



1/2 



1 - 



47r 



1/2 



The free energy is then rewritten as 

A 



£ — Sttk, 



52 



2kc: 



max 



(20) 



(21) 



where a = An/c^g^^ and c = c/c,nax- The free energy 
in Eq. (|21|l has a minimum at c = 0. This minimum is 
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FIG. 11: Net interface length versus time for the case of closed 
vesicles with (pB = 0.5. Top and bottom curves correspond 



to Cj^'"''' and Cj''"'' , respectively. The slope of the dotted line 
is 0.5. 



,(50) 




absolute if the area of the domain is smaller than ag = 
47r(K/A)^. Otherwise, the free energy is lower for c > 
0. These calculations imply that for a given k/X, the 
onset of domain capping occurs when their average radius 
exceeds Rq = 2k/ X, in qualitative agreement with our 
simulation results. 

In order to verify the arguments above, the bending 
modulus and line tension for the cases of a^s = 50 and 
lOOe are extracted numerically. Details of the numerical 
approach used to derive these quantities are described 
in Appendices A and B. We obtain a bending modulus 
K « 8.4e, and a line tension A = 7.1e/rc and 5Ae/rc- The 
number of domains at the onset of capping is calculated 
as 



ttRI 



(22) 



where Av is the vesicles area. In the following table, the 
number of domains at the onset of capping, Nc, obtained 
using Eq. 122|l and that obtained from the simulations 
are shown. 







Ar 


Nc (from Eq. ^) 


91 


54 


Nc (from simulation) 


71 


38 



It is interesting to note that our numerical results are 
fairly in good agreement with the theoretical values. The 
discrepancy is reasonable considering the fact that in 
Eq. H21I) . only lowest order terms, in curvature, are ac- 
counted for. In the simulation, however, a cap is not 
uniformly curved. 

Since domain growth prior to capping is governed by 
a t^^^ law, the onset time of domain capping should be. 




FIG. 12: Snapshot sequences for open membranes with (j)B = 
0.5, with UAB = 50. (I) and (II) correspond to projected 
areas Lx x Ly — 86rc x 86rc and 78rc x 78rc, respectively. 
Snapshots from top to bottom correspond to t — 50r, 200r, 
400r and 800r, respectively. 



tcr ^ {>^/^)^- III the simulation, this time scale is found 
to be approximately 400t and lOOOr for systems Aj^^"-* 

and ^2^'''', respectively. The ratio between these two 
times is very close to the ratio of line tensions given by 

(^A^Gw) / A^aoo) ^ « 0.45. 

In Fig. 8, the mean square of the distance of A-lipid 
heads from the vesicle's center of mass, (Ar)^, is plot- 



ted for systems, A'^'^^\ 



A'^^\ A'^"' , together with the 



^(50) 
.(100) 



case with small excess area AY""". This fi gure clearly 
illustrates that domains capping is shifted towards later 
times as the line tension is decreased. 



C. Late time dynamics of closed vesicles with 
excess area and (j)B = 0.3 

The dynamics in systems ■4^^"°^ and A'^^^'' departs 
from each other after about t — 400r, as shown in Fig. 3: 
The dynamics of domain growth in the case of high excess 
area speeds up at late times, as compared to the case with 
low excess area. As shown in Fig. 3, most domains of 
the minority component in the case with low excess area 
remain flat with a curvature equal to that of the majority 
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FIG. 13: Interfacial length (bold curves) and the square root 
of the second moment (thin curves) versus time for open mem- 
branes with (j>B = 0.5, with uab ~ 50. Red and blue curves 
correspond to a projected area Lx x Ly — 86rc x 86rc and 
78rc X 78rc, respectively. The slope of the dotted line is —1/2. 



average interfacial length, Zc, of a single cap as 




Domain growth proceeds via the Brownian motion of 
domains, leading to their coalescence as they collide. 
Therefore, the mean square of the distance travelled by 
a single domain, 

(f - Dt, (24) 

where D Uc and cP ^ Av/Nc, Ay being the area of 
the vesicle. Using as well, the fact that the net area of 
the _B-domains is a constant of time, i.e., NJ'^ const, 
one then obtains 



and the net interfacial length, Lc — NJ,c is then given by 

- (26) 
in excellent agreement with our numerical results. 



component. Capping, in this case, is suppressed due to 
the lateral tension induced by volume constraint and low 
excess area. In contrast, in the case where excess area is 
high, domains reshape into caps, allowing the interfacial 
length to decrease rapidly. The presence of excess area 
allows some of the caps to further reshape into buds as 
shown in Fig. 2, which then vesiculate, since the line 
tension is relatively high in this case. Bud vesiculation 
occurs over a relatively short period of time, as indicated 
in Fig. 9. We confirm this by performing simulations of 
a vesicle composed of two coexisting domains and with 
high excess area. We found that once the bud is formed, 
it vesiculates within a time period of about IOt. 

The vesiculation of some B-domains results in a 
marked decrease of excess area leading the vesicle with 
the remaining B-domains to acquire a much more spher- 
ical shape. Once all excess area has been released, do- 
main growth crosses over to a regime characterized by 
L{t) ~ with a « 0.4, as shown in Fig. 3. Fig. 4 
demonstrates that the number of domains remaining on 
the vesicle decreases with time as Nb ~ t^^ with a 
growth exponent (3 « 2/3. 

We believe that the higher growth exponent at late 
times, a « 0.4, is the result of coalescence of well formed 
caps following their Brownian motion. Indeed, as Fig. 1 
shows domains of the minority component are more akin 
to hemispherical caps at late times. Assuming that the 
curvature of this caps are determined by the competi- 
tion between the bending energy and the line tension at 
the interface between A and B lipids, we can write the 



D. Case of closed and open membranes with 

4>B = 0.5 

We now focus on the dynamics of phase separation of 
multicomponent closed vesicles with symmetric volume 
fractions of the two components. In Fig. 10, snapshots 
corresponding to the case with high excess area and with 
high line tension are shown. The corresponding inter- 
face length versus time is shown in Fig. 11. This figure 
shows that the characteristic domain size in the case of 
= 0.5 is much more pronounced than that with asym- 
metric composition. Furthermore the growth exponent at 
intermediate times is a « 1/2, larger than that for the 
case of — 0.3. 

We must note that in the case of (j)B = 0.3, where the 
lipid domain structure is circular, two measures where 
used to characterize domain growth. These correspond 
to the interface length and the average cluster size. In 
the case of (j>B = 0.5, only the interface length was so far 
presented. In order to investigate the robustness of the 
corresponding growth exponent, a = 1/2, we also per- 
formed simulations of open membranes extending along 
the xy-plane. These simulations allow us to extract an- 
other length scale from the composition structure factor, 
S{q,t) = {\c^^{t)\^)/L,Ly, as R{t) = 27r/g2, where 



92 (t) 



J dqS{q,t) 



1/2 



(27) 



The effect of tension, and thus the influence of bend- 
ing modes on phase separation, can also be investigated 
through varying the projected area, L^Ly. 
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In Fig. 12, sequences of snapshots of open membranes 
with (jiB = 0-5 and projected areas L^Ly = (80 x 80)r^ 
and (78 x 78)r^ are displayed. Corresponding interface 
length, L{t), and the square root of the second moment, 
92 (i), are presented in Fig. 13. It is clear from this figure 
that both lengths scale as with a « 1/2. 

The exact mechanism leading to this exponent is not 
clear at present. A growth exponent, a = 1, has been 
predicted in the case of two component fluids with inter- 
connected structures in three-dimensions and is the result 
of the instability of the tubular domains against the peri- 
staltic modes [23 ■ Such instability does not exist in pure 
two-dimensional fluids [2^. On the other hand, sev- 
eral simulations on purely two-dimensional fluids based 
on molecular dynamics [23, model H [2^ and lattice 
gas [23, have seen a growth exponent, a = 1/2: A the- 
oretical argument for this exponent is however lacking. 
The lipid bilayer, being a two-dimensional fluid embed- 
ded in a three-dimensional solvent, is clearly more com- 
plicated than a purely two- or three-dimensional fluid. 



IV. SUMMARY AND CONCLUSIONS 

In summary, we presented a detailed study of the phase 
separation dynamics of self-assembled bilayer fluid mem- 
branes, with hydrodynamic effects, using dissipative par- 
ticle dynamics. We considered both open and closed 
membranes, and investigated the effect of composition, 
line tension and surface tension. In all cases, hydrody- 
namics is found to affect the coarsening dynamics at all 
time scales. In the case of closed vesicles with off-critical 
quenches, rich dynamics was observed, with crossovers 
depending strongly on line tension, area-to- volume ratio. 
The early dynamics in this case, is governed by the coales- 
cence of small flat patches. In the presence of excess area, 
the later dynamics dynamics is characterized by the coa- 
lescence of caps. The crossover between the two regimes 
depends strongly on line tension, and includes an inter- 
mediate vesiculation regime for high enough line tension. 
In the case of tensed vesicles, no crossover is observed in 
the dynamics. In the case of critical quenches, the growth 
dynamics is qualitatively different and no crossovers were 
observed. 



Acknowledgments 

The authors would like to thank O.G. Mouritsen and 
L. BagatoUi, M. Rao and G.I. Menon for useful discus- 
sions and critical comments. MEMPHYS is supported by 
the Danish National Research Foundation. Part of the 
simulations were carried out at the Danish Center for 
Scientific Computing. This work is supported in part by 
the Petroleum Research Fund of the American Chemical 
Society. 



Appendix A: Numerical Extraction of the Bending 
Modulus 



The bending modulus is extracted from the power 
spectrum of the long-wavelength out-of-plane fluctua- 
tions of an open lipid membrane, extending along the 
xj/-plane, in a box with periodic boundary conditions 
along the three directions. The Hclfrich Hamiltonian of 
a membrane, in terms of the principal curvatures, ci(r) 
and C2(r), is given by 



H{C1,C2) 



da 



-(ci+C2-2co) +KC1C2 + 



(28) 

where a, k and k are the tension, the bending modulus 
and the Gaussian bending modulus of the membrane re- 
spectively. Co is its spontaneous curvature. In the case 
of a one-component membrane, and if the lipid lateral 
densities of the two leaflets are equal, the spontaneous 
curvature vanishes. Furthermore, if the membrane con- 
serves its topology, the integral of the Gaussian term is 
independent of the membrane conformation. If we as- 
sume that the height of the membrane is represented by a 
single- valued function /i(x), and in the case of small fluc- 
tuations, the Hamiltonian can be rewritten in the Monge 
representation as. 



n{h) 



(29) 



The structure of the membrane can then be inferred from 
the structure factor, defined as the Fourier transform of 
the height-height correlation function. 



5.(q) = -i(i^ n 



dxe 



^Mx) V , (30) 



where q — {qx,Qy)- If the higher order powers of h are 
omitted in the Hamiltonian, and after the equi-partition 
theorem is invoked, one finds the the membrane height 
structure factor is given by 



(31) 



The membrane position, /i(x), is defined as the location 
of the median point of the hydrophobic region. In Fig. 14, 
the structure factor for a one component open membrane 
in a box of dimensions = Ly = 86rc and Lz = 40rc 
is shown. The other parameters are exactly the same as 
those presented in Sec. II. The bending modulus and the 
tension on the membrane are obtained from examining 
the structure factor at small wavevectors. The tension 
on the membrane is found from the intercept of the a{q) 
vs. q^ curve with the y-axis where. 



C7{q) 2kBT/q^Sh{q). 



(32) 



The bending modulus is obtained from the slope of (j{q), 
at small q's, as illustrated in Fig. 14. We obtained a value 
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FIG. 14: a{q) = 2/q^Sh{q) versus obtained from a sim- 
ulation of a one-component open membrane with projected 
area x Ly — 86rc x SGvc- The slope of the curve at small 
wavevectors determines the bending modulus, and the inter- 
cept with the y-axis determines the surface tension of the 
membrane. 



K « 8e, when fc^T — e. The value of k in our model is 
reasonably in agreement with the experimental values for 
lipid bilayers. 



Appendix B: Numerical Extraction of the Line 
Tension 

The in-plane tension of a membrane extending along 
the xy-plane is calculated from the pressure tensor as 



a — Lz 



1 



{Pxx + Pyy) 



(33) 



where the pressure tensor Pnf i is calculated using the 
Irving-Kirkwood formalism 



(34) 



To calculate the line tension, the membrane is pre- 
pared such that it consists of A and B coexisting phases, 
separated by two interfaces which are parallel to the x- 
axis. The tension of the membrane now contains a two- 
dimensional bulk component plus a one-dimensional con- 
tribution due to the interfaces between the A and B com- 
ponents. The line tension is then calculated from the 
difference between the tension of a membrane with two 
interfaces between the A and B components and that of 
a membrane composed with >l-lipids only. 



(35) 



a and a' were calculated on system sizes with dimen- 
sions Lx = Ly = 86rc and = 40rc. We found that 
A « TAe/rc and 5.5e/rc for a/ijha = 50e and lOOe, re- 
spectively. If we assume that the lipid bilayer thickness 
is 4 nm, then our numerically calculated line tensions 
correspond to 2.3 x 10~^''J^m and 1.73 x 10~^''J/im, re- 
spectively. 

Although there are no experimental data for the line 
tension between coexisting lipid phases, a simple es- 
timation can be given by A « ^U/l, where AU = 
{Uaa — Ubb)/2 — Uab ^ lOksT, where Uafs are the var- 
ious lipid pair interactions, and I ~ 0.8 nm is the lateral 
length scale associated with a lipid molecule. One then 
finds A ~ 10~^^ 3 ^m, in agreement with our numerical 
values ;3L |. 
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